{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0d3f0e73",
   "metadata": {},
   "source": [
    "# 机器读心术之文本挖掘与自然语言处理第5课书面作业\n",
    "学号：207402  \n",
    "\n",
    "**书面作业**  \n",
    "1. 在上周作业的基础上，体验HMM第二个基本问题的维特比算法，计算上周幻灯片第34页“掷骰子HMM”中可观测序列为“1635273524”，反求出最大可能的隐藏状态序列，可使用任何编程语言或手工计算。  \n",
    "2. 假设语料库为  \n",
    "研究生物很有意思  \n",
    "他大学时代是研究生物的  \n",
    "生物专业是他的首选目标  \n",
    "他是研究生  \n",
    "\n",
    "（1）试直接用肉眼观察加注分词标注（BMES），然后基于此语料建立相应的HMM，求解出全套模型参数  \n",
    "（2）【可选】如果只有原始语料，没有分词标注，可以建立HMM吗？如果能建立，训练过程是怎样的？  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba9f5665",
   "metadata": {},
   "source": [
    "## 作业1\n",
    "在上周作业的基础上，体验HMM第二个基本问题的维特比算法，计算上周幻灯片第34页“掷骰子HMM”中可观测序列为“1635273524”，反求出最大可能的隐藏状态序列，可使用任何编程语言或手工计算。  \n",
    "**解：**  \n",
    "根据宗庆华书中维特比算法定义如下："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e37b167",
   "metadata": {},
   "source": [
    "维特比变量$\\delta_t(x)$是在时间$t$时，HMM沿着某一条路径到达状态$s_i$，并输出观察序列$O_1O_2...O_t$的最大概率：\n",
    "$$\n",
    "\\delta_t(i)=\\underset{q_1,q_2,\\ldots,q_{t-1}}{\\operatorname{max}}P(q_1,q_2,\\ldots,q_t=s_i,O_1O_2\\ldots O_t|\\mu) \\tag{6.21}\n",
    "$$\n",
    "同时设置另外一个变量$\\psi_t(i)$用于路径记忆，让$\\psi_t(i)$记录该路径上状态$s_i$的前一个（在时间$t-1$）状态。  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f859fac",
   "metadata": {},
   "source": [
    "**步1** 初始化：  \n",
    "$$\n",
    "\\begin{align*}\n",
    "\\delta_1(i)&=\\pi_i b_i(O_1) \\quad,\\quad 1 \\le i \\le N  \\\\\n",
    "\\psi_1&=0\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cebf09af",
   "metadata": {},
   "source": [
    "**步2** 归纳计算：\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\delta_t(j)&=\\underset{1\\le i\\le N}{\\operatorname{max}}[\\delta_{t-1}(i)\\cdot a_{ij}]\\cdot b_j(O_t)\\quad ,\\quad 2\\le t \\le T;1\\le j\\le N \\\\\n",
    "\\psi_t(j)&=\\underset{1\\le i\\le N}{\\operatorname{arg \\, max}}[\\delta_{t-1}(i)\\cdot a_{ij}]\\cdot b_j(O_t)\\quad ,\\quad 2\\le t\\le T;1\\le i\\le N\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edb4d7fa",
   "metadata": {},
   "source": [
    "上面关于$\\psi_t(j)$的公式应该不对，后面不需要乘以$b_j(O_t)$:\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\psi_t(j)&=\\underset{1\\le i\\le N}{\\operatorname{arg \\, max}}[\\delta_{t-1}(i)\\cdot a_{ij}]\\quad ,\\quad 2\\le t\\le T;1\\le i\\le N\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e6d773c",
   "metadata": {},
   "source": [
    "**步3** 终结：\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\hat{Q}_T&=\\underset{1\\le i\\le N}{\\operatorname{arg\\,max}}[\\delta_T(i)] \\\\\n",
    "\\hat{P}(\\hat{Q}_T)&=\\underset{1\\le i\\le N}{\\operatorname{max}}[\\delta_T(i)]\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a588ad0",
   "metadata": {},
   "source": [
    "**步4** 路径（状态序列）回溯：\n",
    "$$\n",
    "\\hat{q}_t=\\psi_{t+1}(\\hat{q}_{t+1})\\quad ,\\quad t=T-1,T-2,\\ldots,1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "95090617",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 状态转换概率矩阵定义如下：\n",
    "a=np.array(\n",
    "    [\n",
    "        [1/3.,1/3.,1/3.],\n",
    "        [1/3.,1/3.,1/3.],\n",
    "        [1/3.,1/3.,1/3.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "# 发射概率矩阵定义如下：\n",
    "b=np.array(\n",
    "    [\n",
    "        [1/6,1/6,1/6,1/6,1/6,1/6,  0,  0],\n",
    "        [1/4,1/4,1/4,1/4,  0,  0,  0,  0],\n",
    "        [1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8]\n",
    "    ]\n",
    ")\n",
    "\n",
    "#初始概率分布定义如下：\n",
    "pi=np.array([1/3,1/3,1/3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "6ae1dc72",
   "metadata": {},
   "outputs": [],
   "source": [
    "def viterbi_proc(a,b,pi,seq,N):\n",
    "    '''\n",
    "    维特比算法\n",
    "    a: 状态转换概率矩阵\n",
    "    b: 发射概率矩阵\n",
    "    pi: 初始概率分布\n",
    "    seq: 观察序列\n",
    "    N: 状态数\n",
    "    '''\n",
    "    T=len(seq)\n",
    "    delta=[]\n",
    "    delta.append(pi*b[:,0])#第1步，初始化\n",
    "    psi=[]\n",
    "    \n",
    "    for t in range(1,T,1):\n",
    "        delta_t=np.zeros(N)\n",
    "        psi_t=np.zeros(N)\n",
    "        for j in range(N):\n",
    "            delta_t_sub_1=delta[-1]\n",
    "            delta_t[j]=np.max(delta_t_sub_1 * a[:,j])*b[j,int(seq[t])]\n",
    "            psi_t[j]=np.argmax(delta_t_sub_1 * a[:,j])\n",
    "        delta.append(delta_t)\n",
    "        psi.append(psi_t)\n",
    "    rpsi=[]\n",
    "    rpsi.append(np.argmax(delta[-1]))\n",
    "    psi.pop()\n",
    "    L=len(psi)\n",
    "    for i in range(L):\n",
    "        t=L-i-1\n",
    "        rpsi.append(int(psi[t][int(rpsi[-1])]))\n",
    "    rpsi.reverse()\n",
    "    return np.max(delta[-1]),rpsi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "14ae0fc8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1.1963376510358935e-12, [1, 2, 1, 0, 1, 2, 1, 0, 0])\n"
     ]
    }
   ],
   "source": [
    "print(viterbi_proc(a,b,pi,'1635273524',3))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12aef861",
   "metadata": {},
   "source": [
    "## 2 作业2\n",
    "\n",
    "假设语料库为：  \n",
    "* 研究生物很有意思  \n",
    "* 他大学时代是研究生物的\n",
    "* 生物专业是他的首选目标\n",
    "* 他是研究生  \n",
    "\n",
    "（1）试直接用肉眼观察加注分词标注（BMES），然后基于此语料建立相应的HMM，求解出全套模型参数  \n",
    "（2）【可选】如果只有原始语料，没有分词标注，可以建立HMM吗？如果能建立，训练过程是怎样的？  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b144d11",
   "metadata": {},
   "source": [
    "**解：**  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eb3f579",
   "metadata": {},
   "source": [
    "(1)  \n",
    "先进行肉眼标注：  \n",
    "* $\\overset{B}{研}\\overset{E}{究}\\overset{B}{生}\\overset{E}{物}\\overset{S}{很}\\overset{S}{有}\\overset{B}{意}\\overset{E}{思}$  \n",
    "* $\\overset{S}{他}\\overset{B}{大}\\overset{E}{学}\\overset{B}{时}\\overset{E}{代}\\overset{S}{是}\\overset{B}{研}\\overset{E}{究}\\overset{B}{生}\\overset{E}{物}\\overset{S}{的}$  \n",
    "* $\\overset{B}{生}\\overset{E}{物}\\overset{B}{专}\\overset{E}{业}\\overset{S}{是}\\overset{S}{他}\\overset{S}{的}\\overset{B}{首}\\overset{E}{选}\\overset{B}{目}\\overset{E}{标}$  \n",
    "* $\\overset{S}{他}\\overset{S}{是}\\overset{B}{研}\\overset{M}{究}\\overset{E}{生}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1cf0d0",
   "metadata": {},
   "source": [
    "在本例中，隐藏序列是BMES，共有4个状态即B、M、E、S。观察序列就是上面的中文四句话。我们要求解全套模型参数，即：$\\pi(i),a_{ij},b_j(k)$，其中$1\\le i,j \\le 4$，这里令$s_1$=B, $s_2$=M,$s_3$=E,$s_4$=S。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ca4b711",
   "metadata": {},
   "source": [
    "这里共有汉字21个，因此$1\\le k \\le 21$。这里令：  \n",
    "$$\n",
    "q_1=研, q_2=究, q_3=生, q_4=物, q_5=很, q_6=有, q_7=意 \\\\\n",
    "q_8=思, q_9=他, q_{10}=大, q_{11}=学, q_{12}=时, q_{13}=代, q_{14}=是 \\\\\n",
    "q_{15}=的, q_{16}=专, q_{17}=业, q_{18}=首, q_{19}=选, q_{20}=目, q_{21}=标 \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bec543e0",
   "metadata": {},
   "source": [
    "根据语料库及肉眼标注情况可得：  \n",
    "$$\n",
    "\\pi=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "\\frac{1}{2} & 0 & 0 & \\frac{1}{2}\\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42569e5c",
   "metadata": {},
   "source": [
    "$$\n",
    "a=\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "0 & \\frac{1}{12} & \\frac{11}{12} & 0\\\\\n",
    "0 & 0 & 1 & 0\\\\\n",
    "\\frac{5}{12} & 0 & 0 & \\frac{4}{12}\\\\\n",
    "\\frac{1}{2} & 0 & 0 & \\frac{2}{5}\\\\\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3a13de6",
   "metadata": {},
   "source": [
    "根据宗庆华书中的向前向后算法：  \n",
    "**步1** 初始化：随机地给参数$\\pi(i),a_{ij},b_j(k)$赋值，使其满足如下条件： "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bdd477e",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\sum_{i=1}^N \\pi(i)&=1 \\\\\n",
    "\\sum_{j=1}^N a_{ij}&=1 \\quad,\\: 1 \\le i \\le N \\\\\n",
    "\\sum_{k=1}^M b_j(k)&=1\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69eebbc1",
   "metadata": {},
   "source": [
    "**步2** EM计算：  \n",
    "**E-步骤** 由模型$\\mu_i$根据公式(6.24)和(6.25)计算期望值$\\epsilon_t(i,j),\\gamma_t(i)$。  \n",
    "**M-步骤** 用E-步骤得到的$\\epsilon_t(i,j),\\gamma_t(i)$，根据公式(6.26)(6.27)(6.28)重新估计参数$\\pi(i),a_{ij},b_j(k)$的值，从而得到模型$\\mu_{i+1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04649c44",
   "metadata": {},
   "source": [
    "**步3** 循环计算：  \n",
    "令i=i+1，重复执行EM计算，直到$\\pi(i),a_{ij},b_j(k)$收敛。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b63f7230",
   "metadata": {},
   "source": [
    "$$\n",
    "\\alpha_t(i)=P(O_1O_2\\cdots O_t|q_t=s_i,\\mu) \\tag{6.12}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7523a76",
   "metadata": {},
   "source": [
    "$$\n",
    "\\beta_t(i)=P(O_{t+1}O_{t+2}\\cdots O_T|q_t=s_i,\\mu) \\tag{6.15}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cb16a25",
   "metadata": {},
   "source": [
    "给定HMM $\\mu$和观察序列$O=O_1O_2\\cdots O_T$，在时间$t$位于状态$s_i$，在时间$t+1$位于状态$s_j$的概率$\\epsilon_t(i,j)$为\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\epsilon_t(i,j)&=P(q_t=s_i,q_{t+1}=s_j|O,\\mu)   \\\\\n",
    "&=\\frac{P(q_t=s_i,q_{t+1}=s_j, O|\\mu)}{P(O|\\mu)} \\\\\n",
    "&=\\frac{\\alpha_t(i)a_{ij}b_j(O_{t+1})\\beta_{t+1}(j)}{\\sum_{i=1}^N\\sum_{j=1}^N \\alpha_t(i)a_{ij}b_j(O_{t+1})\\beta_{t+1}(j)} \\tag{6.24}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc45bc89",
   "metadata": {},
   "source": [
    "给定HMM $\\mu$和观察序列$O=O_1O_2\\cdots O_T$，在时间$t$位于状态$s_i$的概率$\\gamma_t(i)$为\n",
    "$$\n",
    "\\gamma_t(i)=\\sum_{j=1}^N \\epsilon_t(i,j) \\tag{6.25}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "507918a5",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\bar{\\pi}_i&=P(q1=si|O,\\mu)=\\gamma_1(i) \\tag{6.26}\\\\\n",
    "\\bar{a}_{ij}&=\\frac{q_i转移到q_j的次数}{q_i转移到任一态的次数 }  \\\\\n",
    "&=\\frac{\\sum_{t=1}^{T-1}\\epsilon_t(i,j)}{\\sum_{t=1}^{T-1}\\gamma_t(i)} \\tag{6.27}\\\\\n",
    "\\bar{b}_j(k)&=\\frac{Q中从状态q_j发出v_k的次数}{Q中到达q_j的次数 } \\\\\n",
    "&=\\frac{\\sum_{t=1}^{T}\\gamma_t(j)\\delta(O_t,v_k)}{\\sum_{t=1}^{T}\\gamma_t(j)} \\tag{6.28}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d72f1d07",
   "metadata": {},
   "source": [
    "实现源代码如下：  \n",
    "先给参数们一个初始值,这里初始值是按均匀分布随意给的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "be913004",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# 状态转换概率矩阵定义如下：\n",
    "a=np.array(\n",
    "    [\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.],\n",
    "        [1/4.,1/4.,1/4.,1/4.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "# 发射概率矩阵定义如下：\n",
    "b=np.array(\n",
    "    [\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.],\n",
    "        [1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.,1/21.]\n",
    "    ]\n",
    ")\n",
    "\n",
    "#初始概率分布定义如下：\n",
    "pi=np.array([1/4.,1/4.,1/4.,1/4.])\n",
    "\n",
    "#定义观察序列\n",
    "q={'研': 0, '究': 1, '生': 2, '物': 3,'很': 4,'有': 5,'意': 6,\n",
    "   '思': 7, '他': 8, '大': 9, '学':10,'时':11,'代':12,'是':13,\n",
    "   '的':14, '专':15, '业':16, '首':17,'选':18,'目':19,'标':20}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "de31c702",
   "metadata": {},
   "outputs": [],
   "source": [
    "def forward_backward_proc(a,b,pi,q,seq,N,K):\n",
    "    delta=0.05\n",
    "    maxturns=1000\n",
    "    T=len(seq)\n",
    "    a_i=a\n",
    "    b_i=b\n",
    "    pi_i=pi\n",
    "    a_i1=a\n",
    "    b_i1=b\n",
    "    pi_i1=pi\n",
    "    turns=0\n",
    "    while True:\n",
    "        #计算alpha变量\n",
    "        alpha=[]\n",
    "        alpha.append(pi_i*b_i[:,q[seq[0]]])\n",
    "        for t in range(T-1): # 循环T-1次\n",
    "            alpha.append(np.matmul(alpha[-1],a_i)*b_i[:,q[seq[t+1]]])\n",
    "        # print('alpha:',alpha)\n",
    "        #计算beta变量\n",
    "        beta=[]\n",
    "        beta.append(np.ones(N))\n",
    "        for t in range(T-1,0,-1): # 循环T-1次\n",
    "            beta.append(np.matmul(a_i,b_i[:,q[seq[t+1-1]]])*beta[-1])\n",
    "        beta.reverse()\n",
    "        # print('beta:',beta)\n",
    "        #计算epsilon变量\n",
    "        e=[]\n",
    "        for t in range(T-1):#从计算公式上看，只能到T-1，到了T就没有意义了？\n",
    "            e_t=np.zeros((N,N))\n",
    "            denominator=0\n",
    "            numerator=0\n",
    "            for i in range(N):\n",
    "                for j in range(N):\n",
    "                    denominator+=alpha[t][i]*a_i[i][j]*b_i[j][q[seq[t+1]]]*beta[t+1][j]\n",
    "            for i in range(N):\n",
    "                for j in range(N):\n",
    "                    numerator=alpha[t][i]*a_i[i][j]*b_i[j][q[seq[t+1]]]*beta[t+1][j]\n",
    "                    e_t[i][j]=numerator/denominator\n",
    "            e.append(e_t)\n",
    "        #计算gamma变量\n",
    "        gamma=[]\n",
    "        for t in range(T-1):#从计算公式上看，只能到T-1，到了T就没有意义了？\n",
    "            gamma.append(np.sum(e[t],axis=0))\n",
    "        #重新计算pi\n",
    "        pi_i1=gamma[0]\n",
    "        #重新计算a\n",
    "        a_i1=np.zeros((N,N))\n",
    "        denominator_i1=np.zeros(N)\n",
    "        for t in range(T-1):\n",
    "            a_i1+=e[t]\n",
    "            denominator_i1+=gamma[t]\n",
    "        denominator_i1=denominator_i1.reshape(N,1)\n",
    "        a_i1=a_i1/denominator_i1\n",
    "        #重新计算b\n",
    "        b_i1=np.zeros((N,K))\n",
    "        for k in range(K):\n",
    "            for j in range(N):\n",
    "                b_tmp1=0\n",
    "                b_tmp2=0\n",
    "                for t in range(T-1):\n",
    "                    b_tmp1+=gamma[t][j]*(q[seq[t]]==k)\n",
    "                    b_tmp2+=gamma[t][j]\n",
    "                b_i1[j][k]=b_tmp1/b_tmp2\n",
    "        #看pi,a,b是否收敛？\n",
    "        turns+=1\n",
    "        if turns>=maxturns:\n",
    "            print('reach max turns!')\n",
    "            print(np.max(np.abs(a_i1-a_i)),np.max(np.abs(b_i1-b_i)),np.max(np.abs(pi_i1-pi_i)))\n",
    "            break\n",
    "        if np.max(np.abs(a_i1-a_i))<=delta and np.max(np.abs(b_i1-b_i))<=delta and np.max(np.abs(pi_i1-pi_i))<=delta:\n",
    "            break;\n",
    "#         if np.linalg.norm(a_i1-a_i)<=delta and np.linalg.norm(b_i1-b_i)<=delta and np.linalg.norm(pi_i1-pi_i)<=delta:\n",
    "#             break\n",
    "    return a_i1,b_i1,pi_i1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "e2f53611",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha: [array([0.01190476, 0.01190476, 0.01190476, 0.01190476]), array([0.00056689, 0.00056689, 0.00056689, 0.00056689]), array([2.6994925e-05, 2.6994925e-05, 2.6994925e-05, 2.6994925e-05]), array([1.28547262e-06, 1.28547262e-06, 1.28547262e-06, 1.28547262e-06]), array([6.12129818e-08, 6.12129818e-08, 6.12129818e-08, 6.12129818e-08]), array([2.91490389e-09, 2.91490389e-09, 2.91490389e-09, 2.91490389e-09]), array([1.38804947e-10, 1.38804947e-10, 1.38804947e-10, 1.38804947e-10]), array([6.60975939e-12, 6.60975939e-12, 6.60975939e-12, 6.60975939e-12])]\n",
      "beta: [array([5.55219789e-10, 5.55219789e-10, 5.55219789e-10, 5.55219789e-10]), array([1.16596156e-08, 1.16596156e-08, 1.16596156e-08, 1.16596156e-08]), array([2.44851927e-07, 2.44851927e-07, 2.44851927e-07, 2.44851927e-07]), array([5.14189047e-06, 5.14189047e-06, 5.14189047e-06, 5.14189047e-06]), array([0.00010798, 0.00010798, 0.00010798, 0.00010798]), array([0.00226757, 0.00226757, 0.00226757, 0.00226757]), array([0.04761905, 0.04761905, 0.04761905, 0.04761905]), array([1., 1., 1., 1.])]\n",
      "reach max turns!\n",
      "0.0 0.09523809523809526 5.551115123125783e-17\n",
      "[[0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25]]\n",
      "[[0.14285714 0.14285714 0.14285714 0.14285714 0.14285714 0.14285714\n",
      "  0.14285714 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.14285714 0.14285714 0.14285714 0.14285714 0.14285714 0.14285714\n",
      "  0.14285714 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.14285714 0.14285714 0.14285714 0.14285714 0.14285714 0.14285714\n",
      "  0.14285714 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]\n",
      " [0.14285714 0.14285714 0.14285714 0.14285714 0.14285714 0.14285714\n",
      "  0.14285714 0.         0.         0.         0.         0.\n",
      "  0.         0.         0.         0.         0.         0.\n",
      "  0.         0.         0.        ]]\n",
      "[0.25 0.25 0.25 0.25]\n"
     ]
    }
   ],
   "source": [
    "af,bf,pif=forward_backward_proc(a,b,pi,q,'研究生物很有意思',4,21)\n",
    "print(af)\n",
    "print(bf)\n",
    "print(pif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70e12fd7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
